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This  paper  presents  2D  wave-current  interaction  model  for  evaluating  nearly  horizontal  wave-induced  currents  in  the  surf-zone 
and  coastal  waters.  The  hydrodynamic  model  is  the  two-dimensional  depth-averaged  nonlinear  shallow  water  equations  by  us¬ 
ing  an  unstructured  non-staggered  and  multiple-level  quadtree  rectangular  mesh,  this  mesh  information  is  stored  in  simple  data 
structures  and  it  is  easy  to  obtain  a  locally  high  resolution  for  important  region.  The  intercell  fluxes  are  computed  based  on  the 
HLL  (Harten-Lax-van  Leer)  approximate  Riemann  solver  with  shock  capturing  capability  for  computing  the  dry-to-wet  inter¬ 
face  of  coastal  line.  The  effects  of  pressure  and  gravity  are  included  in  source  term  in  the  model,  this  treatment  can  simplify  the 
computation  and  eliminate  numerical  imbalance  between  source  and  flux  terms.  The  wave  model  readily  provides  the  radiation 
stresses  that  represent  the  shortwave-averaged  forces  in  a  water  column  for  the  hydrodynamic  model  and  the  wave  model  takes 
into  account  the  effect  of  wave-induced  nearshore  currents  and  water  level.  The  coupling  model  is  applied  to  verify  different 
experimental  cases  and  real  life  case  of  considering  the  wave-current  interaction.  The  calculated  results  agree  with  analytical 
solution,  experimental  and  field  data.  The  results  show  that  the  modeling  approach  presented  herein  should  be  useful  in  simu¬ 
lating  the  nearshore  processes  in  complicated  natural  coastal  domains. 
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1  Introduction 

Reliable  prediction  of  wave  motion  in  coastal  areas  is  cru¬ 
cial  to  coastal  engineering  applications  associated  with 
nearshore  morphologic  change  and  harbor/inlet  mainte¬ 
nance.  In  some  areas,  however,  ambient  tidal  and  other  cur¬ 
rents  can  be  strong  and  their  effect  on  wave  transformation 
can  be  substantial.  On  entering  shallow  coastal  waters, 
waves  are  modified  by  shoaling,  refraction,  diffraction  and 
wave-current  interaction  before  the  wave  energy  dissipates 
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[1,2].  Radiation  stresses  from  the  excess  momentum  flux  of 
the  waves  produce  mean  water  level  changes  (setting  up  and 
setting  down)  and  longshore  or  cross-shore  currents.  These 
wave-induced  currents  and  water  level  setups  play  a  major 
role  in  sediment  and  contaminant  transport,  pollutant  mix¬ 
ing,  etc.  in  the  nearshore  region.  Very  often  the  flow  pattern 
consists  of  the  combined  interaction  of  tidal  currents  and  a 
wave-induced  current  field.  The  complexity  is  further  en¬ 
hanced  by  the  influence  of  the  flow  field  on  wave  propaga¬ 
tion.  In  order  to  protect  as  well  as  develop  these  coastal  ar¬ 
eas,  knowledge  of  wave-current  interactions  and  wave- 
induced  bottom  shear  stress  is  thus  necessary  to  understand 
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sediment  dynamics,  pollutant  mixing  and  transport  and  bot¬ 
tom  properties. 

In  order  to  solve  the  two-dimensional  shallow  water 
equations,  a  recent  trend  is  the  use  of  finite  volume  method 
(FVM)  [3-9].  This  is  due  to  the  simplicity  of  implementa¬ 
tion,  combined  with  good  flexibility  for  space  discretization. 
On  the  other  hand,  FVM  can  be  interpreted  as  a  fi¬ 
nite-difference  method  applied  to  the  conservative  form  of 
the  mechanical  balances  in  arbitrary  coordinates  and  requires 
less  computational  effort  than  a  finite  element  method  [10, 
11].  To  obtain  an  accurate  resolution  of  discontinuities  of 
the  flow  motion,  some  different  schemes  have  been  devel¬ 
oped  by  researchers.  Because  dam-break  flows  and  storm 
wave  flooding  are  usually  in  the  mixed  flow  regimes  and 
with  discontinuities,  the  often  used  numerical  schemes  are 
the  shock-capturing  schemes,  such  as  approximate  Riemann 
solvers  and  TVD  (Total  Variation  Diminishing)  schemes. 
Recently,  significant  progress  in  the  simulation  of  flood 
wave  propagation  due  to  a  dam-break  event  was  achieved 
by  introducing  the  TVD  schemes  [12-14]  and  approximate 
Riemann  solvers  [15,  16].  Many  of  these  methods  have 
been  verified  and  successfully  applied  to  the  simulation  of 
dam-break  flow  over  the  fixed  beds  and  estuary  tidal  flows 
by  two-dimensional  models  [8,  9,  17,  18].  Because  super¬ 
critical  flows  can  be  seen  due  to  storm  waves  and  surge 
tides  through  breached  barrier  islands  or  coastal  inlets,  hy¬ 
drodynamic  models  should  be  capable  of  simulating  multi¬ 
ple  flow  regimes  such  as  subcritical,  transcritical,  or  super¬ 
critical  flows.  The  wave-induced  nearshore  currents  in  a 
coast  make  flow  pattern  more  complex.  An  exclusive  wave 
model  is  necessary  to  couple  with  the  shallow  water  flow 
model  in  order  to  simulate  both  the  wave  field  and  the 
wave-induced  flow  field.  Most  two-dimensional  depth- 
averaged  models  of  wave-current  interaction  are  based  on 
the  fixed  rectangular  grids  or  uniform  curvilinear  bound- 
ary-fitted  grids  and  therefore  simulate  the  large-scale  flow 
features  on  a  structured  grid  system  [19-22],  but  it  is  not 
very  efficient  computationally  for  the  grid  to  be  of  uni¬ 
formly  high  resolution  throughout  the  entire  numerical  do¬ 
main.  Quadtree  grids  have  gained  increasing  popularity  in 
recent  years  due  to  many  of  their  obvious  advantages:  they 
are  cheap  and  automatic  to  generate,  mesh  information  is 
stored  in  simple  data  structures  and  it  is  easy  to  obtain  a 
locally  high  resolution.  Some  coupling  wave-current  models 
with  quadtree  grid  have  been  developed  for  evaluating 
nearly  horizontal  wave-induced  currents  in  the  surf-zone 
[23,  24], 

In  the  present  study,  an  explicit  2D  shallow  water  hy¬ 
drodynamic  model  is  developed  based  on  a  finite  volume 
method  with  quadtree  mesh.  In  order  to  simulate  wave  field 
and  wave-induced  flow  field  in  a  coast  area,  the  hydrody¬ 
namic  model  is  combined  with  a  steady  state  wave  trans¬ 
formation  and  deformation  model  (CMS-Wave).  The  sys¬ 
tematic  integration  of  the  wave  model  with  the  hydrody¬ 
namic  model  helps  to  develop  a  shock-capturing  coastal 


flow  model  that  can  simulate  wave  transformations  and  de¬ 
formations  together  with  complex  flow  scenarios.  The  pro¬ 
posed  model  is  applicable  to  both  wet-  and  dry-bed  and  ex¬ 
tensively  by  predicting  dam-break  flows  and  wave-current 
interactions  for  coastal  area  with  natural  topography.  The 
results  will  be  presented  in  the  next  sections. 


2  Physical  model 

2.1  Wave  model 


The  wave  model  is  a  two-dimensional  spectral  wave  action 
model  [25].  The  model  formulation  is  based  on  the  para¬ 
bolic  approximation  equation  including  diffraction  terms 
and  energy  dissipation  due  to  wave  breaking  and  bottom 
friction.  The  model  simulates  steady  state  spectral  transfor¬ 
mation  of  directional  random  waves  and  waves  can  propa¬ 
gate  from  the  seaward  boundary  toward  shore  and  vice 
versa.  The  wave  model  takes  into  account  the  effect  of  wave 
breaking,  shoaling,  refraction,  diffraction,  wave-current 
interaction,  bottom  friction  and  the  wave-action  balance 
equation  is  as  follows: 


9(CriV)  |  9(cyv)  |  9 (CgN) 
dx  dy  d0 

K  r  9  f  2ad N)  ccg  9 2N 

2<r[9yy  dy  )  2  9 y 

-  shN  -  Sou, 


where  x,  y  are  the  horizontal  coordinates  in  two  directions; 
N=E(a ,  6  )/<r  is  the  wave-action  density  to  be  solved  and  is 
a  function  of  frequency  a  and  direction  9,  which  is  defined 
as  angle  of  wave  relative  to  the  x-axis.  The  spectral  wave 
density  E(a,  9)  represents  the  wave  energy  per  unit  water 
surface  area  per  frequency  interval  in  eq.  (1).  The  first  term 
on  the  right-hand  side  of  eq.  (1)  represents  wave  diffraction 
as  formulated  from  the  parabolic  wave  approximation  as¬ 
sumption.  A  default  value  of  k  =2.5  is  used  for  the  diffrac¬ 
tion  intensity  parameter  suggested  by  Mase  in  the  present 
study  [25].  C  and  Cg  are  the  wave  celerity  and  group  veloc¬ 
ity,  respectively,  and  is  parameterization  of  wave  break¬ 
ing  energy  dissipation.  Sou  is  the  source  terms  (e.g.,  wind 
forcing,  bottom  friction  loss,  nonlinear  wave-wave  interac¬ 
tion  term  etc.)  [26].  These  characteristic  wave  velocities 
with  respect  to  x,  y  and  9  coordinates  are  accordingly  C„  Cy 
and  Cg  defined  as 


Cx  =  Cg  cos  9  +  u,  Cy  =  Cg  cos  9  +  v,  (2) 
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where  u  and  v  are  current  velocity  components  in  the  x  and 
y  directions,  respectively,  k  is  the  wave  number,  and  h  is 
water  depth.  Eq.  (1)  is  solved  for  wave-action  density  N  by 
a  forward-marching  first-order  upwind  finite-difference 
method  with  rectangle  grid,  the  detailed  description  of  solu¬ 
tion  methodology  can  be  found  in  ref.  [25].  The  importance 
of  this  function  is  examined  for  four  wave  breaking  formu¬ 
las,  in  this  study,  the  extended  Miche  formula  is  used  for 
calculating  wave  breaking  term,  a  summary  description  of 
the  parameterized  wave  breaking  formulas  is  presented  by 
Zheng  et  al.  [27]. 


2.2  Wave-induced  nearshore  current  equation 

The  numerical  model  used  consists  of  the  two-dimensional 
shallow  water  equations  describing  the  conservation  of 
mass  and  momentum.  The  shallow  water  equation  written  in 
conservation  and  vector  form  are 


dU  9  E  9  G  9  E,  9  G,  a 

d  '  d-  +  S. 


-  H - h  - 

9 1  dx  9  y 


-  +  - 
dx  9  y 


(4) 


In  which.  U  is  the  vectors  of  conserved  variables,  E(U), 
G(U),  Ed(U)  and  GA(U)  in  eq.  (4)  are  convection  fluxes  and 
diffusion  fluxes  in  the  x  and  y  directions,  respectively,  S(  LJ) 
is  source  terms  which  can  be  defined  respectively  as  follows: 


h 

hu 

hv 

dull 

u  = 

hu 

,  E  = 

huu 

,  G  = 

huv 

>E*  = 

V‘~dT 

hv 

huv 

hvv 

dvh 

dx 


0 

0 

Gd  = 

dull 

dy 

,  5  = 

ghp-Sb,-TSl-Tia  +  fchv 

OX 

1 

CD  1  9^ 

^  I  s- 
1 _ 

Sh^-Sby-^y-T^-fJiu 

(5) 


where  g  is  the  acceleration  due  to  gravity,  .ST,  and  Sby  are 
bed  shear  stress  terms  with  x  and  y  components  defined  by 
the  velocities,  Z  is  water  level,  Zb,  the  bed  surface  elevation 
above  a  reference  datum.  The  unit-widths  hu  and  hv  are  the 
conservative  dependent  variables,  grouped  in  the  column 
vector  U.  rsx  and  rsv  are  wave  force  terms  with  x  and  y  com¬ 
ponents;  fc  is  the  Coriolis  parameter.  rWJ  and  rwv  are  the  sur¬ 
face  wind  stresses,  respectively. 

In  the  depth-averaged  parabolic  model  the  eddy  viscosity 
is  given  for  pure  current: 

v,  =  auji.  (6) 


In  which  u*  is  the  bed  shear  velocity  u,  =  [cf(u 2  +  v2)]1/2 

and  a  is  an  empirical  coefficient  between  0. 3-1.0.  The 
transition  between  surf-zone  mixing  and  oceanic  mixing 
seaward  of  the  breakpoint  is  represented  in  CMS-M2D  [20] 
by  a  weighted  mixing  coefficient  specified  as 


v,  =  (1  -  D)au,h  +  DAumH,  (7) 


where  A  is  the  mixing  parameter  equal  to  1 .0,  um  is  the  am¬ 
plitude  of  the  horizontal  component  of  the  wave  orbital  ve¬ 
locity  at  the  bottom  given  by 


u  = 


gHT 


2 A  cosh 


2n  {h  +  ij'j 

~~X 


(8) 


where  H  is  wave  height,  T  is  wave  period,  >]  is  deviation  of 
the  water-surface  elevation  from  the  still-water  level,  A  is 
the  coefficient  for  um  which  can  be  seen  in  ref.  [20], 
D  =  H/(h  +  rj). 

Sbx  and  Sby  are  bed  shear  stress  terms  with  x  and  y  com¬ 
ponents  defined  by  the  velocities  for  the  situation  without 
waves, 

Sbr  =  gii2u\ju2  +  v2 h  Sbv  =  gn2v\lu2  +  v2 h  ^ ,  (9) 


where  n  is  the  manning’s  roughness  coefficient.  In  the  case 
of  simulating  wave-induced  currents,  the  near  bed  orbital 
velocity  has  to  be  taken  into  account  here.  Nishimura’s  ap¬ 
proximation  is  used  to  calculate  the  shortwave-averaged 
bottom  stresses  under  combined  currents  and  waves  in  this 
study  [28]; 


S=n2hA< 


Sby=n2h~'4< 


a>„ 


t/„c+— cos  y 


u  + 


'  EL 


cos  ^  sin  y 


r 


cos/ sin  ^ 


u  + 


V  , 


CO,  r\ 

U  + — —cos’  v 
U 

V  UW  c  J 


(10) 


v  . 


ctHs 

01  = - £ - . 

k  sinh  (kH) 


(11) 


In  which  y  is  wave  angle  relative  to  the  x-axis,  Hs  is  signifi¬ 
cant  wave  height.  [/wc  is  the  resultant  velocity  from  the 
shortwave-averaged  and  the  wave-orbital  velocities  given 
by 


Uwc  =  i  Jtr  +  v2  +  col  +  2  («  cos  y  +  v  sin  y)  cob  | 

—^u2  +  v2  +  cob  —  2 (u cosy  +  vsin y^ cob  .  (12) 


+  — , 
2 


The  wave  force  is  calculated  by  using  the  following  equa¬ 
tions; 


dx  9  y 


=  ~ 


dx  9  y 


(13) 


(14) 


where  Sxx,  Sxy,  and  Syy  are  wave-driven  radiation  stresses. 
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Radiation-stress  tensor  calculations  are  based  on  linear 
wave  theory.  They  represent  the  summation  of  standard 
tensor  formulations  across  the  defined  spectrum.  For  a  co¬ 
ordinate  system  with  the  x-axis  oriented  normal  to  the 
shoreline,  the  tensor  components  are  as  follows  [20]: 


S„  = 


r[E(a,0) 

f 

0.5 

JJ  2 

V 

1+  2  k(h  +  n)  icos220  +  l)-O.5 

sivA\2k  [h  +  rj  )  y  ' 


AoAO, 

(15) 


S„  = 


// 


E{a,9) 


0.5 


n  2k{h  +  rj) 
sinh  2k  (h  +  ?]j 


sin  20 


AaAO,  (16) 


5W  = 


rr  ,  J  (  2  k(h  +  rj)},  , 

E(a,0 )  0.5  1  + -  ' \  (sin2  20  + l)-0. 

JJ  v  ’  y  sinh  2k  (h  +  rj)  Jv  7 

The  Coriolis  parameter  is  given  as 
fc  =  2Qsin(^), 


dcrdt?. 

(17) 

(18) 


where  f2  is  the  angular  frequency  of  the  Earth’s  rotation, 
and  (j)  is  latitude. 

Surface  wind  stresses  are  given  by  the  following  empiri¬ 
cal  wind  shear  stress  formulations: 

C,,  =cd— W2sin(/?),  (19) 

Av 


Evy  =ci— W2  cos (P), 

Pw 


(20) 


where  Cd  is  the  wind  drag  coefficient;  pa  and  pw  are  density 
of  air  and  water;  W  is  wind  speed  measured  at  10  m  above 
the  sea  level;  and  /?  is  the  angle  of  wind  direction  relative  to 
x-axis. 


has  one  or  two  faces  on  each  of  its  north,  south,  west,  and 
east  sides,  as  shown  in  Figures  1  and  2,  so  that  the  computa¬ 
tional  mesh  will  be  less  complicated  [29]. 

The  data  structure  for  the  quadtree  mesh  can  be  managed 
in  several  ways:  block-structured,  hierarchical  tree,  and 
fully  unstructured.  The  block-structured  approach  divides 
the  domain  into  multiple  blocks,  each  of  which  is  treated  as 
structured.  Interfaces  between  blocks  need  to  be  specially 
handled  to  ensure  mass  and  momentum  balance  through 
them.  The  tree  data  structure  uses  parent  and  child  relations 
and  requires  tree  traverse  to  determine  the  mesh  connec¬ 
tivity.  In  the  fully  unstructured  approach,  all  the  cells  are 
numbered  in  a  one-dimensional  sequence,  and  pointers  are 
used  to  determine  the  connectivity  of  neighboring  cells  for 
each  cell.  Among  the  three  approaches,  the  fully  unstruc¬ 
tured  approach  is  simpler  and  thus  is  used  in  this  study.  As 
mentioned  in  Section  1,  another  issue  in  simulation  of  in¬ 
compressible  flow  is  the  location  of  primary  variables:  ve¬ 
locity  and  water  level.  On  the  non-staggered  grid,  all  the 
primary  variables  are  located  at  the  center  of  cells.  The 
computer  code  based  on  non-staggered  grid  is  simpler  and 
requires  less  memory,  and  the  non-staggered  grid  is  simpler 
in  handling  the  interface  between  coarse  and  fine  cells. 


3  Grid  type  and  numerical  discretization 

3.1  Quadtree  grid  and  data  structure 

Because  of  the  complexity  of  computational  domain,  a  sim¬ 
ple  structured  rectangular  mesh  requires  a  large  number  of 
cells  to  resolve  the  detailed  flow  pattern  near  the  dam 
breaking  location,  navigation  channels  and  in-stream  struc¬ 
tures.  To  optimize  the  use  of  computational  resources,  we 
use  the  multiple-level  quadtree  rectangular  mesh  with  local 
refinement.  In  this  mesh,  various  levels  of  fine  cells  are 
placed  close  to  the  dam-break  locations  and  in-stream 
structures  or  inlet  areas  where  the  flow  gradients  are  high, 
while  coarse  grids  are  used  in  the  low-gradient  regions.  To 
simplify  the  mesh,  a  cell  is  refined  by  splitting  into  four 
equal  child  cells.  Corresponding  to  this  refining,  any  cell 
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Figure  2  Control  volume  in  a  quadtree  mesh  based  on  non-staggered 
grid. 
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3.2  Finite  volume  method 


The  numerical  representation  of  physical  domain  is  ob¬ 
tained  by  using  an  unstructured  mesh,  composed  of  quad¬ 
tree  elements.  Each  of  them  must  be  considered  as  an  ele¬ 
mentary  control  volume.  Eq.  (4)  can  be  integrated  over  a 
control  volume  V  as 


f  !E+f  ^dv+f^dv-f 

J  v  dt  J  v  3 x  J  v  3v  Ji 


3  E 


d  G 


dy 


9 A 

h  dx 


dy 


f  3  Gb  r 

— ^dy+  SdV.  (21) 

J  V  3v  J  v 


<v  dy 


Applying  the  Gauss  theorem  to  eq.  (21),  discretizing  eq.  (21) 
over  a  time  interval  by  using  the  explicit  scheme  in  time, 
one  can  obtain: 
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The  detailed  definition  about  notations  (e,  w,  n  and  s)  can  be 
found  in  Figure  2,  the  interpolation  for  water  level  (velocity 
etc.)  at  between  cell  (node  P)  and  a  child  cell  (node  El)  can 
be  expressed  in  eqs.  (23)  and  (24): 


Jaxp2+A y2 

fin  -  V  I — ==r=. 

(23) 

Ay£,  +y]Ax;+A yEi 

Zei  =  fin*  zEi  +(1  -  fin)*  zp. 

(24) 

3.3  HLL  scheme 


A  detailed  description  of  the  HLL  (Harten-Lax-van  Leer) 
scheme  can  be  found  in  Toro  [17]  and  Ying  [8]  including  a 
complete  discussion  regarding  the  Riemann  problem  and 
the  reasons  for  a  special  treatment  of  the  cells  located  on  the 
wave  front  or  shore  line  (wet/dry  boundary).  The  HLL 
scheme  assumes  and  defines  the  flux  at  an  interface  as 
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In  which  f)  and  FR  are  the  flux  evaluated  at  the  left-hand 
and  the  right-hand  sides  of  each  cell  interface.  F  denotes 
the  flux  at  the  intermediate  state,  given  by 


F  = 


SrF  l  A^r  +  AlAr  (Gr  Gl) 


Ar  Al 


(26) 


In  which  t/L  and  UR  are  the  conservative  variable  vectors 
evaluated  on  the  left-hand  and  the  right-hand  sides  of  each 
cell.  The  symbols  Al  and  Sr  represent  the  celerity  of  the 
waves,  separating  constant  states  of  the  local  Riemann 
problem  solution  at  cell  interfaces,  they  can  be  estimated  as 
follows: 


SL  =  min -  yfgh[, u  - yfgli^j,  (27) 

SR  =  min  (mr  +  yfgh^,  u  +  yjgh~  ) ,  (28) 


The  symbols  hL  and  /;R  are  the  water  depth  of  the  left  and 
right  states.  Note  that  for  a  dry  bed  problem,  the  wave 
speeds  Sl  and  SR  can  be  estimated  according  to  the  follow¬ 
ing  expressions: 

SL=uL-  SghR ,  SR=uL  +  2 y]ghL  for  right  dry,  (31) 
SL  =  uR  -  2^J ghR ,  SR  =  uR  +  yj ghR  for  left  dry  bed.  (32) 

3.4  Stability  analysis 

A  variable  Ar,  adapted  to  hydraulic  parameters  variability,  is 
used.  It  is  well  known  that  computation  stability  for  explicit 
scheme  model  requires  a  Courant-Friedrichs-Lewy  (CFL) 
number  less  than  1 .  This  parameter  in  the  FVM  context  de¬ 
fined  as 

CFL  =  Max  +  yfgh  )>^“(|v|  +  v/s^)|  <  1-  (33) 

The  time  step  can  be  either  automatically  selected  by  the 
program  based  on  the  consideration  of  the  CFL  condition  or 
directly  specified  by  the  user. 


4  Model  validation  and  verification 

4.1  Dam  break  on  dry  and  wet  beds 

The  aim  of  this  test  case  is  to  study  the  ability  of  the  model 
to  simulate  discontinuous  flow  and  the  front  wave  propaga¬ 
tion  over  dry  and  wet  beds.  Considering  a  horizontal  and 
frictionless  channel  of  1200  m  long  and  1  m  wide,  a  column 
of  water  is  separated  by  a  dam  located  at  a- 500  m.  This  test 
case  is  simulated  for  both  wet  and  dry  bed  situations,  the 
initial  upstream  water  depth  is  10  m.  The  downstream  water 
depth  is  5  m  for  wet  bed  and  0  m  for  dry  bed,  the  tolerance 
depth  for  dry  bed  is  set  as  0.0001  m.  The  space  discretiza¬ 
tion  is  set  as  Ax=  1  m  and  time  step  is  set  as  At-  0. 1  s,  at  t= 
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0  the  dam  is  removed  instantaneously.  The  results  are  re¬ 
ported  up  to  30  s  after  the  removal  of  the  dam.  Figures  3(a) 
and  3(b)  compare  the  water  depth  and  the  velocity  for  the 
wet  bed  case  and  Figures  4(a)  and  4(b)  compare  water  depth 
and  the  velocity  for  the  dry  bed  case.  The  computed  water 
depth  and  velocities  follow  the  exact  solutions. 


0  300  600  900  1200 

Distance  (m) 

Figure  3  Dam  break  on  wet  bed.  (a)  Water  depth;  (b)  velocity. 
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Figure  4  Dam  break  on  dry  bed.  (a)  Water  depth;  (b)  velocity. 


4.2  Flow  over  a  triangular  obstacle 

This  test  is  carried  out  to  reproduce  the  laboratory  dam- 
break  flow  over  a  triangular  hump  recommended  by  the  EU 
CAD  AM  project  [15].  The  physical  experiment  included 
complex  hydraulic  properties  such  as  shocks,  transitions 
between  wet  and  dry  beds  and  flow  over  an  obstacle.  The 
laboratory  setup  is  illustrated  in  Figure  5.  The  experiment 
consists  of  a  reservoir  with  water  up  to  0.75  m  contained  by 
a  dam  at  x  =15.5  m  and  a  dry  bed  downstream  within  a  rec¬ 
tangular  channel  of  20.5  m  in  length.  A  symmetric  triangu¬ 
lar  obstacle  (6  m  long,  0.4  m  high)  is  placed  on  the  channel 
with  its  peak  located  at  13  m  downstream  of  the  dam.  In 
order  to  observe  depth  evolutions  gauge  points  are  located 
at  4  m  (G4),  10  m  (G10),  11  m(Gll),  13m(G13)and  20  m 
(G20)  from  the  dam  as  shown  in  Figure  5.  The  upstream 
boundary  is  wall  and  the  downstream  boundary  is  free  flow. 
The  Manning’s  roughness  coefficient  is  set  as  0.0125  s  m_1/3 
with  the  same  value  of  the  experiment.  In  the  computation, 
v=0.1  m,  t  =0.01  s  and  /?=0.0001  m  as  dry  bed  are  used  and 
the  simulation  is  carried  out  for  100  s.  The  water  depth 
variations  with  respect  to  time  at  the  five  gauge  points  are 
shown  in  Figures  6(a)-6(e).  It  can  be  seen  that  the  predicted 
water  depth  evolutions  and  arrival  time  of  the  wave  are 
quite  comparable  with  the  measured  data  at  every  gauge 
point.  The  transition  from  wet  to  dry  at  the  gauge  point  G13 
is  well  predicted.  The  general  trend  is  captured  well  at  G20 
(after  the  hump)  and  the  water  depth  is  slightly  overesti¬ 
mated  before  20  s,  which  is  also  predicted  by  other  re¬ 
searchers  by  using  different  numerical  schemes  [14,  15]. 
Overall,  the  comparison  between  the  numerical  predictions 
and  measurements  is  satisfactory,  the  wet-dry  transitions  are 
simulated  with  very  high  accuracy  by  the  present  model. 
This  confirms  the  effectiveness  of  this  method  on  flux  and 
the  effects  of  pressure  and  gravity  in  source  term.  All  the 
results  are  computed  by  using  the  depth-averaged  2D  model, 
the  reason  with  slight  discrepancies  may  be  that  the  2D 
models  cannot  accurately  simulate  the  complex  3D  flow 
phenomenon. 

4.3  Wave-induced  longshore  current 

Kuriyama  and  Ozaki  [30]  carried  out  field  measurements  of 
the  longshore  current  at  Hazaki  Oceanographical  Research 
Facility  (HORF)  located  on  the  Japan  Pacific  coast.  The 
HORF  research  pier  is  427  m  long,  and  the  current  meas¬ 
urements  were  made  from  the  HORF  pier  by  using  a  float. 
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Figure  5  Geometry  and  gauge  locations  in  the  experiment. 
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Figure  6  Computed  and  measured  water  depth  variations  at  stations  (a)  G4;  (b)  G10;  (c)  Gil;  (d)  G13;  and  (e)  G20. 


The  beach  at  HORF,  having  a  mean  slope  of  1/60,  often 
includes  several  longshore  bars,  leading  to  complex  wave 
transformation  with  shoaling,  breaking,  and  reforming  tak¬ 
ing  place  as  shown  in  Figure  7.  The  significant  wave  height 
is  2.6  m  and  significant  wave  period  is  8.86  s,  the  wave  di¬ 
rection  angle  a  is  27°  in  this  case  study  [1],  A  time  step  of 
0.2  s  is  used,  the  Manning’s  n  is  set  as  0.015  and  the  simu¬ 
lation  stops  until  steady  state  reaches.  Figure  8  shows  the 
calculated  and  measured  distribution  of  the  longshore  cur¬ 
rent.  Breaking  on  the  seaward  side  of  the  bar  yields  two 
peaks  in  the  longshore  current  distribution,  which  is  in 
agreement  with  the  measurements.  The  seaward-most  peak 
agrees  with  the  observed  one,  including  the  correct  location, 
the  shoreward-most  peak  has  the  correct  magnitude,  but  is 
located  somewhat  shoreward  of  the  measured  peak.  The 
mixing  parameter  is  set  to  A=  1 .0  in  this  study.  Figure  9 
compares  the  calculated  and  measured  significant  wave 
height.  Figure  10  shows  the  calculated  longshore  current 
pattern  for  Kuriyama  and  Ozaki  field  experiment  [30].  The 
results  demonstrate  that  the  coupling  model  gives  predic¬ 
tions  in  close  agreement  with  experimental  data,  whereas 
numerical  model  somewhat  overestimates  the  wave  heights 
in  the  breaker  region.  As  would  be  expected,  the  wave 


Figure  7  Bed  profile  for  Kuriyama  and  Ozaki  field  experiment  [30]. 


Figure  8  Calculated  and  measured  longshore  current  for  Kuriyama  and 
Ozaki  field  experiment  [30]. 
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Figure  9  Calculated  and  measured  wave  height  for  Kuriyama  and  Ozaki 
field  experiment  [30]. 


Figure  10  Calculated  longshore  current  for  Kuriyama  and  Ozaki  field 
experiment  [30]. 


heights  rapidly  decrease  to  zero  onshore  after  breaking.  This 
rapid  loss  in  height  of  the  broken  waves  may  promote 
strong  radiation  stresses,  which  in  turn  may  create  large  and 
possibly  unstable  currents. 

5  Model  applications 

5.1  Tidal  flow  through  an  idealized  inlet 

The  developed  model  is  tested  in  a  hypothetical  case  of  tidal 
flow  through  an  idealized  inlet.  The  inlet  is  150  m  wide.  Its 
east  end  is  connected  with  a  rectangular  bay  of  about  1250  m 
x  2750  m  in  x  and  y  directions,  and  its  west  end  is  open  sea. 
The  sea  bed  has  a  constant  equilibrium  cross-shore  profile 
with  a  maximum  still  water  depth  of  about  20  m,  while  the 
bay  and  the  inlet  have  a  flat  bed  that  is  3  m  below  the  still 
sea  level.  Two  jetties  are  located  on  the  sea  side  of  the  inlet. 
A  four-level  quadtree  rectangular  mesh  (dots  are  cell  centers) 
is  used,  as  shown  in  Figure  1 1 .  The  finest  grid  spacing  near 
the  inlet  is  12.5  m  by  12.5  m,  while  the  coarsest  grid  spac¬ 
ing  is  100  m  by  100  m  in  the  sea  and  50  m  by  50  m  inside 
the  bay.  An  M2  tide  is  specified  at  the  sea  and  two 
cross-shore  boundaries,  the  computational  time  step  is  2  s. 
Figure  12  shows  the  simulated  flow  pattern  during  flood 
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Figure  11  Mesh  near  an  idealized  inlet  (dots  are  cell  centers). 


Figure  12  Simulated  flow  pattern  in  the  bay  during  flood  tide  (a)  and  ebb 
tide  (b). 
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and  ebb  tides.  One  can  see  that  the  flood  jet  eddies  are 
simulated  reasonably  well.  The  recirculation  flows  cover 
four  levels  of  the  refined  mesh,  and  the  transition  between 
flows  on  fine  and  coarse  grids  is  very  smooth,  the  flow  pat¬ 
terns  are  symmetric  about  the  centerline  of  the  inlet  and  two 
symmetric  vortexes  can  be  found  beside  jetties.  This  test 
demonstrates  that  the  numerical  discretization  and  solution 
methods  in  the  developed  model  are  adequate  to  handle  the 
unstructured  quadtree  mesh.  Because  there  are  no  meas¬ 
urement  data,  the  performance  of  the  developed  model  is 
only  qualitatively  verified  in  this  hypothetical  case.  Quanti¬ 
tative  validation  of  the  model  in  case  of  unsteady  flow  will 
be  carried  out  in  the  next  case. 

5.2  Tidal  flow  in  Gironde  Estuary 

The  Gironde  Estuary  is  located  in  southwestern  France.  It 
receives  runoff  from  the  Garonne  River  and  the  Dordogne 
River,  and  empties  into  the  Atlantic  Ocean,  as  shown  in 
Figure  13.  The  water-surface  width  varies  from  2  to  14  km, 
and  the  flow  depth  in  the  navigation  channel  is  about  6-30 
m.  The  estuary  is  partially  mixed  and  macrotidal,  with  a  12 
hour  and  25  minutes  tidal  lunar  period  and  a  tidal  amplitude 


of  1.5-5  m  at  the  mouth  [31].  The  simulation  domain  is 
about  80  km  long,  starting  from  the  mouth  to  the  Garonne 
River  and  the  Dordogne  River.  Because  the  domain  is  sim¬ 
ple,  a  uniform  mesh  is  used  here,  with  a  size  of  250  m  x  125 
m  for  each  cell.  The  measured  data  from  May  19  to  25  of 
1975  are  used  to  validate  the  model.  The  computational 
time  step  is  5  s.  At  the  estuary  mouth,  the  tidal  elevation  is 
given  according  to  the  recorded  time  series  at  the  station 
“Pointe  de  Grave”.  At  the  two  upstream  ends,  the  flow  dis¬ 
charges  of  the  Garonne  River  and  the  Dordogne  River  are 
specified  according  to  the  measured  data  at  La  Reole  and 
Pessac.  The  Manning’s  n  is  set  as  0.015.  Figure  14  shows 
the  comparison  of  the  measured  and  simulated  water  levels 
at  four  selected  stations.  The  amplitude  and  phase  are  well 
predicted  by  the  present  numerical  model,  no  obvious  phase 
difference  exists  between  the  measured  and  simulated  tidal 
levels.  Figure  15  shows  the  comparison  of  the  measured  and 
simulated  flow  velocities  at  five  selected  stations.  The 
measured  flow  velocities  are  1  m  under  the  water  surface  and 
1  m  above  the  seabed,  respectively.  The  simulated  depth- 
averaged  flow  velocity  stays  between  them,  the  agreement  is 
reasonably  good.  The  flow  fields  in  flood  and  ebb  tides  are 
reasonably  well  predicted,  as  shown  in  Figure  16. 


Figure  13  Sketch  of  Gironde  Estuary,  France. 
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Figure  14  Computed  and  measured  water  surface  level  variations  at  different  stations. 
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Figure  15  Computed  and  measured  velocities  at  different  stations. 
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5.3  Tidal  flow  in  Grays  Harbor 

Grays  Harbor  is  located  on  the  southwest  Washington  coast 
about  72  km  north  of  the  Columbia  River.  The  estuary  is 
one  of  the  largest  in  the  continental  United  States.  As  part  of 
the  U.  S.  Army  Corps  of  Engineers  Grays  Harbor  Estuary 
Physical  Dynamics  Study,  data  of  current  and  waves  were 
measured  during  September  to  November  of  1999  [1,  32]. 
The  estuary  has  a  wetted  surface  area  of  approximately  146 
km2  at  mean  higher  high  water  and  45  km2  at  mean  lower 
low  water.  The  main  input  of  fresh  water  is  from  the  Cheha- 
lis  River.  The  5  km  wide  entrance  has  two  convergent  rock 
jetties  that  extend  from  spit  points,  as  shown  in  Figure  17. 
The  spectral  waves  from  the  NOAA  buoy  46029  were  input 
data  at  the  model  boundaries  every  3  h.  The  mesh  is  refined 
around  the  jetties  and  near  the  channels,  the  finest  grid 
spacing  is  25  m  near  the  jetties  and  the  coarsest  one  is  800 
m  near  the  offshore  boundary  at  deep  water.  The  computa¬ 
tional  time  step  is  1  s.  The  measured  water  levels  from  the 
station  nearest  to  the  offshore  boundary  are  used  as  the 
boundary  condition.  The  wave  radiation  stresses  calculated 
by  the  wave  model  are  considered  in  the  simulation  of  tidal 
current.  Figure  18  compares  the  computed  and  measured 
water  levels  at  Stations  1,  5  and  tide  2,  and  Figure  19  com¬ 
pares  the  computed  and  measured  current  speeds  at  Stations 
2  and  5  for  a  period  of  8  days  in  the  late  September  of  1999. 
Figure  20  compares  the  computed  and  measured  current 
direction  at  Station  5  with  good  agreement  and  no  obvious 
phase  difference.  Figure  21  compares  the  computed  and 


Figure  16  Simulated  flow  patterns  in  Gironde  Estuary. 


Figure  17  Topography  and  measurement  stations  at  Grays  Harbor. 


Figure  18  Measured  and  calculated  tide  levels  at  Grays  Harbor. 


measured  wave  heights  at  the  selected  stations.  From  these 
figures,  one  can  see  that  the  agreement  between  those  cal¬ 
culated  results  and  measured  data  is  generally  good.  This 
model  can  also  simulate  the  wetting  and  drying  processes 
on  the  floodplain  due  to  tide  change.  This  demonstrates  that 
the  coupling  model  can  handle  the  moving  coastal  line  effi¬ 
ciently  and  reproduce  the  tidal  water  level  and  current  speed 
reasonably  well. 
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Figure  20  Measured  and  calculated  current  direction  at  Grays  Harbor. 


6  Conclusions 

A  depth-averaged  two-dimensional  tide  and  wave  coupling 
model  is  presented  for  simulating  flow  phenomena  of 
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Figure  21  Measured  and  calculated  wave  heights  at  Grays  Harbor. 


coastal  waters  in  this  study.  Major  characteristics  and  new 
idea  of  this  coupling  model  include: 

1)  This  hydrodynamic  model  applies  the  non-staggered 
grid,  all  the  primary  variables  are  located  at  the  center  of 
cells,  so  the  computer  code  is  simpler  and  requires  less 
memory,  and  the  non-staggered  grid  is  simpler  in  handling 
the  interface  between  coarse  and  fine  cells. 
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2)  Such  an  unstructured  model  with  quadtree  mesh  al¬ 
lows  refining  the  grid  for  important  computational  domain, 
such  as  coastal  inlet  location,  navigation  channel. 

3)  The  intercell  fluxes  are  evaluated  based  on  the  HLL 
approximate  Riemann  solver  with  shock  capturing  property 
for  capturing  accurately  the  moving  waterline  and  therefore, 
capable  of  simulating  multiple  flow  regimes,  in  which  sub- 
critical,  transcritical  or  supercritical  flows  may  happen  in 
the  coastal  waters. 

4)  The  spectral  wave  action  model  can  simulate  wave 
fields  by  accounting  for  wave  breaking,  shoaling,  refraction, 
diffraction,  wind  effect  and  current  effect  in  coastal  zones. 
Four  parameterized  wave  breaking  formulations  are  in¬ 
cluded  in  order  to  calculate  accurately  energy  dissipation 
rate  in  the  wave  model  under  different  conditions. 

5)  Dynamic  coupling  of  the  two  models  requires  ex¬ 
change  of  the  information  in  an  iterative  process,  the  wave 
model  needs  to  provide  the  results  of  radiation  stresses, 
wave  height  and  wave  period  for  the  hydrodynamic  model. 
Meanwhile,  the  hydrodynamic  model  needs  to  provide  the 
water  level  and  velocity  for  the  wave  model. 

The  calculated  results  show  good  performance  of  the 
coupling  model  in  predicting  dam  break  flows,  the  unsteady 
tidal  flow  in  estuaries  and  the  longshore  current  generated 
by  waves.  These  results  show  that  the  developed  model  is 
expected  to  be  a  useful  tool  for  the  study  of  many  coastal 
phenomena  in  the  future. 
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